Ab-initio phase diagram of ultracold 87 Rb in an one-dimensional two-color superlattice 
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We investigate the ab-initio phase diagram of ultracold 87 Rb atoms in an one-dimensional two-color super- 
lattice. Using single-particle band structure calculations we map the experimental setup onto the parameters 
of the Bose-Hubbard model. This ab-initio ansatz allows us to express the phase diagrams in terms of the ex- 
perimental control parameters, i.e., the intensities of the lasers that form the optical superlattice. In order to 
solve the many-body problem for experimental system sizes we adopt the density-matrix renormalization-group 
algorithm. A detailed study of convergence and finite-size effects for all observables is presented. Our results 
show that all relevant quantum phases, i.e., superfluid, Mott-insulator, and quasi Bose-glass, can be accessed 
through intensity variation of the lasers alone. However, it turns out that the phase diagram is strongly affected 
by the longitudinal trapping potential. 

PACS numbers: 67.85.Hj; 03.75.Lm; 67.85.-d 
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I. INTRODUCTION 



Ultracold atomic gases in optical lattices have been a topic 
of active research for about a decade now. One of the research 
thrusts is the use of these systems as experimental quantum 
simulators for a variety of lattice models and allow for de- 
tailed investigations of strongly correlated quantum systems 
in a perfectly controllable environment iH 01 • By tuning the 
laser intensities of the optical lattice alone, one can seamlessly 
drive a system through quantum phase transitions like the su- 
perfluid to Mott-insulator transition O, |4|] . In so-called two- 
color superlattices, additional lasers are used to introduce ir- 
regular lattice topologies which give rise to exotic quantum 
phases like the Bose-glass phase IS Si- 

Strongly correlated particles in periodic potentials are well 
described by Hubbard-type models. Together with powerful 
many-body methods this allows for theoretical studies of the 
phase diagram of ultracold atomic gases in optical lattices ||7|- 
12411 . However, a one-to-one comparison between experiment 
and theory has rarely been done so far because these theoret- 
ical studies usually adopt the generic parameters of the Hub- 
bard model to span the phase diagram. Such a phase diagram 
of ultracold bosonic atoms in a two-color superlattice is shown 
inFig.HIa). 

In this work we establish a closer link to experiments by 
computing the phase diagrams with respect to the natural 
experimental control parameters, which are the intensities 
(S2, s\) of the two lasers generating the one-dimensional opti- 
cal superlattice. Such an experiment- specific phase diagram is 
shown in Fig. QIb). In order to predict this type of phase dia- 
gram we start with single-particle band structure calculations 
to extract the Hubbard parameters for a specific experimen- 
tal setup. Then, the many-body problem is solved using the 
density-matrix renormalization-group (DMRG) algorithm. In 
the following, our band structure plus DMRG approach is in- 
troduced and benchmarked. We discuss the phase diagram of 
a specific experimental setup motivated by Refs. ||5|,|6|] with a 
focus on its dependence on the transverse trapping frequency 
oj ± and the longitudinal trapping frequency co x . 



ID BOSE-HUBBARD MODEL AND BAND STRUCTURE 
CALCULATIONS 



II. 



The single-band Bose-Hubbard model [25] is a widely used 
framework for studying the ultra-low temperature physics of 
strongly correlated, neutral atoms in sufficiently deep optical 
lattices. We assume a one-dimensional lattice with / sites and 
N bosonic atoms. For each site we define the creation (anni- 
hilation) operators aj (0.) with respect to the localized Wan- 
nier states corresponding to the lowest Bloch band. The mean 
occupation-number at each lattice site is given by hi 
The Bose-Hubbard Hamiltonian 



a a.. 



h = Yj {- j um ( aj + i h i + h \ h M ) 

+ -U i (h i -l)h i + e i h\ (1) 



accounts for three basic processes: the tunneling of atoms to 
adjacent sites, the on-site two-body interaction, and the on- 
site potential energy. The site-dependent Hubbard parameters 
Ji,i+i* Uu and 6; define the relative strengths of the individ- 
ual terms and contain all information about depth and topol- 
ogy of the optical potential, and the interaction between the 




1.8 


(b) 








/ // [ 


1.6 










/ / I :; " 


1.4 










v^i3 


1.2 










^-^r^ 


Z 1 










{BG)~s 


0.8 










a^ 


0.6 










^—-^ 


0.4 
0.2 




(SF) 




t 


?t 


..«« 




I 4 


6 


8 


10 


12 14 16 



FIG. 1: (color online) Contour plots of the energy gap computed 
using DMRG for a commensurate superlattice with / = N = 30 as 
(a) a function of the generic Hubbard parameters and (b) a function 
of the experimental laser intensities. The labels mark the domains 
of the superfluid (SF) phase, the homogeneous Mott-insulator (MI) 
phase, and the quasi Bose-glass (BG) phase (taken from Ref. [22]). 



atoms. The phase diagrams spanned directly by these param- 
eters, typically using J iti+ i = /, JJ[ = U and some ansatz for 
€t to account for superlattice structures, are extensively dis- 
cussed in Refs. Il3-[l6i [T9M2TI1 . More recently, the on-site 
energies 3 were calculated directly from the parameters of the 
optical superlattice to provide a closer connection to experi- 
ment I2I . 

In this work, our aim is a discussion of the phase diagram 
using the experimental parameters directly and not the generic 
Hubbard parameters. To this end, an explicit treatment of the 
underlying single-particle physics is necessary. Therefore, we 
start from the optical potential generated by two orthogonal 
polarized standing-wave laser-fields with wavelengths X\ and 
/I2 and the respective potential depths s\ and s^. Furthermore, 
we consider an additional harmonic potential with frequency 
(jj x accounting for the intensity variation of the optical lattice 
through the focusing of the laser beams and a magnetic trap- 
ping potential. Using the recoil energy E r . = ^2 of atoms 

with mass m as a natural energy scale and a phase shift be- 
tween the standing waves, the potential along the x-axes reads: 



V(x) = s\E n sin 



2tt 

— x + 1 
A\ 



+ s 2 E r2 sin z I — x 



-\--moj 2 x 2 . 



(2) 



Throughout this work we consider a setup defined by A2 = 800 
nm and S2 for the primary laser generating the optical lattice 
potential and A\ - 1000 nm and s\ for the secondary laser 
generating the two-color superlattice topology with a phase 
shift of cp = 7r/4.This leads to the commensurate superlattice 
that was also used in previous publications lll41l5[ll9U22ll . 

Before we are able to extract the Hubbard parameters for a 
given potential, we have to determine the localized Wannier 
functions via a single-particle band structure calculation. For 
a periodic potential (s\ - 0, co x = 0) with / sites, we numeri- 
cally obtain the solutions for the Bloch functions i//k(x) in the 
lowest energy band. The quasimomenta are quantized with 
respect to the size / of the optical lattice and are labeled by k 
with k = 0,1,- — ,1-1. A Fourier transformation of the Bloch 
functions with respect to the quasimomenta in the subspace of 
the lowest energy band leads to the Wannier functions 



Wi(x) 



±Y* k (x)e-*»&. 



V/^ 



(3) 



k=0 



The arbitrary phases (fk are chosen such that the resulting 
Wannier functions are maximally localized at their individ- 
ual lattice site i. Using these maximally localized Wannier 
functions, the Hubbard parameters are obtained via the matrix 
elements of the individual terms of the real-space Hamiltonian 



-Ju = 



r * t n 2 d 2 \ 

J^ w ' (x H"^^ +v(x) r (x) 



(4) 
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TABLE I: Calculations of the Hubbard parameters U and J and the 
energy gap AS between first and second Bloch band for a homoge- 
neous lattice (s\ = 0, cj x = 0) with A 2 = 800 nm, a> ± = 2n x 17 kHz, 
a s = 109 r Bo hr, and mass m of 87 Rb . 



The contact interaction is defined by the three-dimensional s- 
wave scattering length a s . The transverse directions are in- 
tegrated out assuming Gaussian wavefunctions with frequen- 
cies coy = qj z = qj_ l |2|]. This one-dimensional description 
is valid as long as the tunneling in the transverse directions 
is strongly suppressed, i.e., as long as the laser intensities in 
these directions are sufficiently large. Following Refs. i^,M 
we consider a gas of 87 Rb atoms with s-wave scattering length 
a s = 109 rBohr and we assume a transverse trapping frequency 
cl)j_ = 2nx 17 kHz. In the last part of this work we will discuss 
the changes in the phase diagram induced by a different value 

Of (JL> ± . 

As a fir/st application of our band structure calculations we 
validate the single-band approximation in the Hubbard model. 
As a by-product from the calculations of the Hubbard parame- 
ters we obtain from the single-particle band structure calcula- 
tion the energy gap A8 between the first and the second band 
of Bloch functions. In Table [J we list some values of AS to- 
gether with the Hubbard parameters U and / in the relevant 
parameter range. Since AS is always about one order of mag- 
nitude larger than U and /, we conclude that excitations to 
energetically higher Bloch bands induced by tunneling or in- 
teraction can be neglected even for shallow optical lattices. 

As a second application of our band structure calcula- 
tions we check for the validity of the restriction to nearest- 
neighbor tunneling and on-site two-body interactions. For the 
weakly and the strongly interacting regime, we calculate the 
respective matrix elements of the Hubbard Hamiltonian using 
Eqs. © and Utj - 2co 1 fia s J dx \Wi(x)\ 2 \Wj(x)\ 2 for the in- 
teraction term. The results are shown in Table HH Even in the 
weakly interacting regime (S2 = 2), the nearest-neighbor tun- 
neling exceeds more-distant tunneling processes by at least 
one order of magnitude. The interaction matrix element for 
neighboring lattice sites is already two orders of magnitude 
smaller than the on-site interaction matrix element. In the 
strongly interacting regime (S2 = 10) we already have two 
orders of magnitude between J^+x and Jtj+2 and five orders 
of magnitude between U(j and Uij+\. Since we focus on the 
intermediate and strong interaction regime, the restriction to 
Jij+i and Ut is well justified. 

So far we have discussed the limit of a homogeneous op- 
tical lattice. As soon as the secondary laser which gener- 
ates the superlattice, or an additional harmonic potential are 
taken into account, a straight-forward band structure calcu- 
lation is no longer doable, because Bloch functions are only 
defined for strictly periodic potentials. Therefore, in order to 
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TABLE II: Higher order tunneling and interaction energies for a ho- 
mogeneous lattice (s\ = 0, cd x - 0) with A 2 = 800 nm, oj ± -2nx\l 
kHz, a s = 109 r Bo hr, and mass m of 87 Rb . 



extract site-dependent Hubbard parameters also for an inho- 
mogeneous lattice we are limited to an approximate scheme 
to obtain localized Wannier functions. We use two different 
approaches to extract the site-dependent Hubbard parameters. 

As a simple ansatz, we consider the secondary laser as a 
perturbation of the strong primary laser (s\ <*c s^). The Wan- 
nier functions are extracted from a conventional band struc- 
ture calculation for a homogeneous lattice defined by the pri- 
mary laser alone. In this approximation the Wannier functions 
are identical for each lattice site. Using these Wannier func- 
tions the Hubbard parameters of each site of the superlattice 
are computed. The site-dependence of the parameters thus re- 
sults exclusively from the superlattice potential V(x) entering 
into the matrix elements © and not from a site-dependence 
of the Wannier functions themselves. As a result, the param- 
eter Ut characterizing the on-site interaction remains constant 
for all lattice sites. An exemplary set of site-dependent Hub- 
bard parameters calculated in this scheme is shown in Fig. [2] 
Please note that we always subtract a global energy constant 
from the Hamiltonian to set 6^ = min{£;} = 0. 

In a more sophisticated scheme we determine the site- 
dependent Wannier functions individually for each site of the 
inhomogeneous lattice using a standard band structure calcu- 
lation for a periodic lattice with a lattice amplitude defined 
by the local depth of the inhomogeneous potential at that par- 
ticular site. In this way, the shape of the Wannier functions 
depends nontrivially on the local structure of the superlattice 
potential. The only reason why the set of Wannier functions 
determined in this way cannot be considered as an exact set 
of localized basis functions results from the minimal viola- 
tion of orthogonality for the Wannier functions of neighbor- 
ing sites. Their mutual overlap is nonzero but always below 
1% in the parameter regime considered in all our calculations. 
Using these individual localized Wannier functions all site- 
dependent Hubbard parameters are computed without further 
approximations. An exemplary set of results is also shown in 

Fig. El 

The comparison of the site-dependent Hubbard parameters 
€t and Jij+i obtained by the two schemes shows very little dif- 
ference. This leads to the conclusion that the second scheme 
provides a sufficiently accurate description of the Hubbard pa- 
rameters in the parameter range under consideration, simply 
because the change induced by first, much cruder approxima- 
tion is small. 




FIG. 2: Site-dependent Hubbard parameters obtained from band 
structure calculations for a two-color superlattice. Simple ansatz 
(gray symbols) and calculations with individual Wannier functions 
(black symbols), both for s 2 = 10 and s\ - 1. Lines to guide the eye. 



The dominant effect on the Hubbard parameters induced by 
the superlattice is the spatial variation of the on-site energies 
€f. This is in agreement with the approximation of the super- 
lattice through this parameter alone Ill4fl^fl9l42lll . However, 
also the tunneling matrix element Jij+i, essentially depending 
on the height of the potential barrier between sites i and i + 1 , 
varies significantly. The on-site interaction matrix element Ut 
exhibits only a weak variation which is introduced by the site- 
dependence of the Wannier functions in our second scheme. 



III. DENSITY-MATRIX RENORMALIZATION-GROUP 

We solve the many-body problem associated with the Bose- 
Hubbard Hamiltonian via the density-matrix renormalization- 
group (DMRG) algorithm I 27l 12811 which is among the most 
powerful quasi-exact methods available for one-dimensional 
lattice models. The so-called infinite- size DMRG algorithm 
is based on an iterative growing procedure. The algorithm is 
schematically depicted in Fig. [3j The individual steps are: 
(0) We start with a block composed of 4 sites and up to 
TVb particles described in a Fock space 9% of dimension D\>. 
Since the Hubbard Hamiltonian conserves the particle num- 
ber, the matrix representation of the block Hamiltonian has a 
block-diagonal form. Each block of the matrix corresponds 
to a Hilbert space with a fixed particle number. (1) To the 
block we attach an additional lattice site with up to N s par- 
ticles described in a Fock space T^ to build the Fock space 
of the system Ts ys = 3% ® T$ with dimension Z) sys = J\D % . 
Again the matrix representation of the system Hamiltonian is 
block diagonal. (2+3) In order to simulate a larger lattice, the 
system is coupled to an analogously constructed environment 
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^en V of dimension 
Auper which is projected to a fixed total particle number, sat- 
isfying N/I = 1 in our case. (3) The ground state | i//q) is 
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FIG. 3: Sketch of the DMRG cycle. The black dot in step 5 marks 
the additional site that was attached without increasing the dimension 
of the Hilbert space of the block. For details see text. 



obtained by diagonalizing the superblock Hamiltonian where 
one can exploit the sparseness of the Hamilton matrix and use 
efficient Lanzcos or Jacobi-Davidson algorithms. (4) The re- 
duced density-matrix is formed by tracing out the environment 
p red = Tr env | <Ao)(<Ao I • (4+5) The Z4 eigenvectors of the re- 
duced density-matrix for the largest eigenvalues are used to 
span the Fock space for a new block Th of length 4 = 4 + 1- 
These eigenvectors build a non-unitary transformation matrix 
O which is employed to construct the new block Hamiltonian 
74 = O^H sys O. All operators coupling the system to the en- 
vironment (which will later couple the new block to the new 
site) and all observables have to be transformed accordingly. 
This cycle is repeated until the final length of the lattice is 
reached. 

The key feature of this algorithm lies in the use of the eigen- 
vectors with largest eigenvalues of the reduced density-matrix 
as a new, truncated basis for the new block. One can show 
that this procedure yields an optimized wavefunction, gives 
the best approximation to expectation values of observables, 
and preserves a maximum of entanglement between system 
end environment IJ28I1 . 

The error in the DMRG algorithm is due to the loss of in- 
formation during the non-unitary basis transformations. It can 
be estimated by summing up the eigenvalues of the discarded 
eigenvectors. A smaller sum consequently means a smaller 
loss of information. In addition one has to consider the re- 
striction to a maximal number of particles max{fzj per lattice 
site. In the complete Hilbert space this would be equal the 
total number of particles N. In general, the stronger the corre- 
lations between the particles, i.e., the larger U, the smaller the 
sum of the residual eigenvalues and the better the approxima- 
tion. 

If disorder is introduced, then only at the very last step of 
the growing procedure the full information about the super- 



lattice topology is available to the Hamiltonian. This leads 
to a poor approximation of the ground state when using the 
infinite- size algorithm only. As an improvement the finite- size 
DMRG is applied. After a complete run of the infinite-size al- 
gorithm up to the desired length of the lattice, the length of 
the superblock is kept fixed and the system grows on the ex- 
pense of the environment and vice versa. During a back and 
forth sweeping, the superlattice topology is sampled while the 
Hamiltonian always takes the whole lattice into account. The 
sweeping continues until all observables are converged. 

The way the transformation matrices O are constructed is 
not uniquely defined by the DMRG algorithm. We would like 
to emphasize that the reduced density-matrix p red is block di- 
agonal and each block has a well defined particle number. One 
can either use the eigenvectors with the largest eigenvalues 
for each subspace of p red , or one can strictly use the first Z4 
eigenvectors with the largest eigenvalues not accounting for 
the block-diagonal structure of p red . In the first scheme one 
might discard eigenvectors with sizable eigenvalues if the re- 
spective subspace has reached its preassigned dimension. In 
the second scheme one might discard a complete subspace of a 
certain particle number in case it has no eigenvector with cor- 
responding eigenvalue among the largest Z4 eigenvalues. If in 
a subsequent step of the finite size algorithm this subspace be- 
comes important again, this might prevent the algorithm from 
converging to the proper ground state. This can be overcome 
by adding noise to the transformation matrices O during the 
first few sweeps in the finite- size algorithm with the aim of 
recovering lost subspaces again [29] . As a third strategy one 
can keep at least one or a few eigenstates from each subspace 
even if their eigenvalues are not among the largest Z4 eigen- 
values. We employ the first strategy because it is technically 
very convenient. However, we checked individual eigenspec- 
tra of p red and conformed that none of the discarded eigenvec- 
tors had sizable eigenvalues. 



IV. OBSERVABLES 

In this section, we introduce the set of observables we em- 
ploy to distinguish the different quantum phases. 

Maximum Number Fluctuation. The number fluctuation 
at lattice site i is given by the variance of the occupation num- 
ber 



«)-<*) 2 



(5) 



The number fluctuation provides information about the local 
mobility of the atoms in the optical lattice. In order to re- 
duce the amount of information we only consider the maxi- 
mum number fluctuation through the lattice 



c^max = maxfcr,-; 



(6) 



Condensate Fraction. In order to determine the fraction 
of atoms that undergo Bose-Einstein condensation we adopt 
the Onsager-Penrose criterion [30] and calculate the natural 
orbitals via the eigensystem of the one-body density-matrix 



Pi] ~ W^j)' The l ar S est eigenvalue N c of the one-body 
density-matrix is associated with the number of condensed 
atoms and defines the condensate fraction 



Jc ~ N 



(7) 



Visibility. In the experiment, most information about the 
atoms in the optical lattice is extracted from the interfer- 
ence pattern obtained by the time-of-flight method. The in- 
terference pattern 1(6) is intimately connected to the quasi- 
momentum structure of the many-body state and can be calcu- 
lated from the Fourier transformation of the one-body density- 
matrix id 



I(S) 






Vij ' 



(8) 



The visibility of the interference fringes v is obtained from the 
maxima and minima of the interference pattern 



y = 



max{ 1(6)} - min{ 1(6)} 
max{ 1(6)} + min{ 1(6)} 



(9) 



Energy Gap. Measuring the excitation spectrum of the 
system also provides a sensitive probe for the different quan- 
tum phases. In the experiment one employs two-photon Bragg 
spectroscopy via an intensity modulation of the optical lat- 
tice. The width of the central interference peak is used as 
a measure of the energy transfer into the atomic cloud 0]. 
The detailed structure of the excitation spectrum has been in- 
vestigated theoretically Ill7l - l20ll . Basic information about the 
excitation spectrum is given by the energy gap AE, which is 
the minimum amount of energy needed to excite the system. 
It is defined by the difference between the energy of the first 
excited state and the ground state 



A£ = £i 



(10) 



An inherent complication in the DMRG framework is the 
calculation of observables. This is because a DMRG calcula- 
tion does not yield eigenstates of the Hamiltonian in a simple 
occupation-number basis representation which could be used 
to compute observables directly. Rather, the matrix represen- 
tations of all observables have to be dragged through all the 
cycles of the DRMG algorithm, i.e., they have to undergo all 
the lossy non-unitary basis rotations. This is why we thor- 
oughly test our DMRG calculations for convergence. 



V. BENCHMARK OF THE DMRG ALGORITHM 

A. Convergence 

Before we employ the DMRG algorithm to compute phase 
diagrams for realistic lattice sizes and particle numbers, we 
have to assess the precision of the numerical DMRG results. 
We follow a twofold strategy. 
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TABLE III: Different bases used for studying the convergence of the 
DRMG calculations. See text for details. 



First, we compare results for the various observables ob- 
tained by DMRG calculations with results from an exact di- 
agonalization scheme (Tj, LL5|] for a small system with / = 
N = 10, where the latter calculations are feasible. In the 
phase diagram shown in Fig.QJa) we observed an error of the 
DMRG calculation below 1% for all observables at U/J > 3 
already for a small DMRG basis with dimension Auper = 338. 
The complete Hilbert space used in the exact diagonalization 
scheme has a dimension of D = 92378 for the / = N = 10 
system. 

Second, in order to validate the results of our DMRG cal- 
culations for larger lattices, where no exact calculations in the 
complete Hilbert space are available, we study the dependence 
of the DMRG results on basis sizes and particle number trun- 
cations used in the algorithm Q2VL Um • If the results for all 
observables do not change while the bases size is increased 
further, the calculation is converged to the exact result. The 
different basis sets we employ are summarized in Table [III 
where max{nj is the maximum number of particles per lattice 
site included in the basis. 

For all following calculations we applied three sweeps in 
the finite-size algorithm. For simplicity we consider straight 
lines through the parameter plane shown in Fig. UX&)- 

Superfluid to Mott-insulator (e max = 0). The upper row 
of images in Fig. |4] shows the observables across the super- 
fluid to Mott-insulator phase transition calculated using the 
three bases specified in Tab. Hm Since the DMRG algorithm 
is tailored to describe strongly correlated systems, we expect 
better agreement of the three different calculations with in- 
creasing U/J. Apart from the energy gap this is confirmed by 
our calculations. Only for U/J < 3 we observe small differ- 
ences between the calculations for <x max and f c . For all values 
of U/J the energy gap is slightly larger when employing the 
DMRG-A basis. This is because we do not explicitly target 
at the first excited state for the calculation of the energy gap. 
Although the ground state has already converged even for the 
small DMRG-A basis, the first excited state needs a larger ba- 
sis to converge as well. 

Mott-insulator to quasi Bose-glass (U/J = 30). We al- 
ready pointed out the importance to use the finite- size DMRG 
algorithm in order to obtain a converged ground state espe- 
cially when irregularities in the optical lattice are considered. 
The results of the observables through the Mott-insulator to 
quasi Bose-glass transition depicted in the lower panel of 
Fig. [4] show that the finite- size algorithm is perfectly con- 
verged for all values of e max /J already for the DMRG-A basis. 

Since, in this manuscript the focus is on the regime of inter- 
mediate and strong interactions, we conclude from our find- 
ings that already the DMRG-A basis is suitable to approxi- 
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FIG. 4: Benchmark of the convergence of the DMRG calculations. From left to right: energy gap AE, maximum number fluctuation cr m2LX , 
condensate fraction f c , and visibility v for a I = N = 30 lattice. Upper panels: SF-MI transition at fixed e max /J = 0. Lower panels: MI-BG 
transition at fixed U/J - 30. All plots show three lines corresponding to a different basis choice: DMRG-A (dotted), DMRG-B (dashed), and 
DMRG-C (solid). 



mate all observables with sufficient precision. Nevertheless, 
we decided to use the larger DMRG-B basis for all following 
calculations. Calculations based on this basis are still numer- 
ically feasible on a desktop PC while providing good results 
also in the weakly interacting limit. 



B. Finite-Size Analysis 

We also have to address the dependence of the observables 
observables on the size of the system. Current experiments 
typically have between 1.5 • 10 4 to 2 • 10 5 atoms in the optical 
lattices QlaL i.e., roughly between 25-60 atoms in each one- 
dimensional array of sites. Thus, we will compare DMRG 
calculations for / = TV = 10, / = N = 30, and / = N = 60 all 
using the DMRG-B basis. In analogy to the above discussion, 
we show plots through the superfluid to Mott-insulator transi- 
tion and the Mott-insulator to quasi Bose-glass transition. 

a. Superfluid to Mott-insulator (e max /J = 0). The results 
are shown in the upper row of Fig.[5j By definition, the maxi- 
mum number fluctuation is a local observable which is calcu- 
lated at one individual lattice site and is, therefore, practically 
independent of the size of the lattice. The energy gap as well 
as the visibility show only small differences between the small 
and the two larger lattices indicating a minor dependence on 
length of the lattice for those observables. However, the con- 
densate fraction depends systematically on the size of the lat- 
tice. The larger the lattice is, the steeper is the decrease of f c 
around U/J ~ 5. One can easily show that for U/J^oo and 
I = N the condensate fraction scales like f c oc \/l 111211 which 
is in-line with our calculations. 



We also performed an additional calculation for the large 
I = N = 50 lattice using the DMRG-C basis. These results 
are not shown in the plots because there are no sizable devi- 
ations to calculations with the DMRG-B basis. Only for the 
condensate fraction at U/J < 3 the DMRG-C basis yields 
slightly larger values, e.g. f c = 0.68 instead f c = 0.63 at 
U/J= 1. This indicates the slower convergence of the DMRG 
algorithm in the weakly interacting regime. The results for all 
other observables remain completely unchanged when going 
to the larger DMRG-C basis. 

b. Mott-insulator to quasi Bose-glass (U/J = 30). The 
lower row of Fig. [5] reveals that the energy gap as well as the 
maximum number fluctuation do not change with the size of 
the lattice across the Mott-insulator to quasi Bose-glass tran- 
sition. The condensate fraction exhibits the previously men- 
tioned 1/7 scaling which is characteristic for large values of 
U/J. The visibility in the small lattice is again slightly smaller 
compared to the two larger lattices. 

Considering this analysis, we conclude that calculations in- 
cluding N = 30 particles on I = 30 lattice sites are sufficient 
to describe realistic experiments. Firstly, because this system 
size is right in the experimental range. And secondly, when 
going to larger systems, there are only small and predictable 
changes for the condensate fraction whereas all other observ- 
ables remain unchanged. 



VI. AB-INITIO PHASE DIAGRAMS 

After the validation of our framework we now discuss the 
experiment- specific phase diagram of an ultracold 87 Rb gas 




\£imLx/J = 




"0 10 20 30 40 50 60 

£max/J 



10 20 30 40 50 60 

£max/J 



1 
0.8 
0.6 
~0.4 
0.2 





5 10 15 20 25 30 
U/J 

U/J = 30 



^ 



1 


"^\ ^max/J = 


0.8 


\ 


0.6 


V 


0.4 


\^^ 


0.2 
n 





10 20 30 40 50 60 

Emax/J 



1 

0.8 
0.6 
0.4 
0.2 





5 10 15 20 25 30 

U/J 

U/J = 30 




10 20 30 40 50 60 

£max/ J 
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with scattering length a s = 109 rBohr in an optical lattice with 
wavelength /fe = 800 nm. The superlattice topology is gen- 
erated by an additional laser with wavelength X\ = 1000 nm 
and relative phase shift of = n/4. The respective optical 
potential depth resulting from the two lasers are given by the 
dimensionless parameters S2 and s\ . The remaining transverse 
lasers of the optical trap enter via the transverse trapping fre- 
quency o) ± which is chosen to be 2n x 17 kHz. Initially, the 
longitudinal trapping frequency co x is set to Hz. 

We have already used these parameters in Fig. [1] to compare 
the experiment- specific phase diagram spanned by S2 and s\ 
with a generic phase diagram spanned by U/J and e max /J ne- 
glecting the site dependence of U and /. Both panels of Fig. [T] 
show the energy gap AE for / = N = 30 obtained from a 
DMRG calculation using the DMRG-B basis. Since the vari- 
ation of ^2 and s\ affects all Hubbard parameters simultane- 
ously, the (s2, s\) phase diagram is distorted in comparison to 
the (U/J, 6 max //) phase diagram. However, the (s^, s\) phase 
diagram reveals that all relevant quantum phases are accessi- 
ble through the variation of the intensity of the two longitudi- 
nal lasers alone, while keeping the other parameters fixed. 

A detailed analysis of the phase diagram for this set of pa- 
rameters is given in Figs.[6ta)-(d), where we show the energy 
gap, the condensate fraction, the maximum number fluctua- 
tion, and the visibility, respectively. 

The superfluid (SF) phase is characterized by a vanishing 
energy gap, large condensate fraction, large number fluctua- 
tions, and maximum visibility. Although we do not compute 
the most stringent order parameter for the SF phase — the 
superfluid fraction Jsl [lj, |23|] — the above signatures allow 
us to identify the SF phase in the region of small S2 up to 
^2 < 6 in the whole range of s\ shown here. Due to the shal- 



low optical potential in this region the tunneling term in the 
Hubbard Hamiltonian (Q]) dominates. This results in a coher- 
ent many-body state which is a prerequisite for the SF phase. 
For S2 = 6 along < s\ < 2 the mean interaction energy is 
U/J ~ 4.5 which explains the presence of the SF phase in the 
whole range of s\. 

In a homogeneous lattice (s\ = or e max /J = 0) a transition 
from the SF phase to the homogeneous Mott-insulating (MI) 
phase occurs around U/J^5 (Tj, |23|] which corresponds to 
^2 = 6.25. This is in-line with our results, because around 
S2 ~ 6 the energy gap steeply increases while the condensate 
fraction, the number fluctuations, and the visibility decrease. 
At ^2 = 16 and si = the ratio of U/J is 60 and the system is 
deep in the homogeneous MI phase showing the characteristic 
large energy gap and vanishing number fluctuations, conden- 
sate fraction, and visibility. 

If we now increase s\ at fixed ^2 = 16, the modulation of 
the site-dependent Hubbard parameters grows rapidly and at 
s\ « 0.6 the spread of the on-site energies becomes compa- 
rable to the average interaction energy, i.e., e max /J « U/J. 
Thus, despite the strong repulsive interaction, it becomes ad- 
vantageous to move an atom from a site with large on-site 
energy to an already occupied site with small on-site energy. 
Due to this redistribution of particles the homogeneous MI 
phase is broken up and the transition to the quasi Bose-glass 
(BG) phase occurs. The commensurate superlattice defined 
by X<i - 800 nm, X\ - 1000 nm, and = n/4 exhibits only 
5 different on-site energies. This small set of on-site energies 
leads to extended domains in the phase diagram. Two of these 
domains are visible in Fig.^a). Only in the transition region 
between them the energy gap vanishes. 

The genuine Bose-glass phase occurs only in an infinite lat- 
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FIG. 6: (color online) Phase diagram in terms of energy gap AE/J, condensate fraction f c , maximum number fluctuation <r max , and visibility 

y for / = N = 30, a s = 109 r Boh r, u± = 2tt x 17 kHz, and cj x = Hz. 
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TABLE IV: Comparison of the Hubbard parameters to analyze effect 
of the longitudinal trapping potential. The parameters are: A 2 = 800 
nm, A\ = 1000 nm, = tt/4, s\ = 0, a) ± = 2n x 17 kHz, and mass 
and scattering length of 87 Rb. 



tice with random on-site energies. It is marked by a com- 
pletely vanishing energy gap. Intuitively this results from a 
continuous distribution of on-site energies permitting the con- 
struction of excited states by moving particles to sites with in- 
finitesimally larger on-site energies associated with infinites- 
imally small excitation energies. We have approached this 
limit using an incommensurate superlattice in a previous pub- 
lication 1221. 



VII. LONGITUDINAL TRAPPING FREQUENCY co x 

Since the aim of this manuscript is the calculation of an 
experiment specific phase diagram for a realistic experimen- 



tal setup, it is compulsory to consider an additional mag- 
netic trapping potential and the intensity variation of the op- 
tical lattice through the focusing of the laser beams. To this 
end we have introduced a harmonic potential with frequency 
(jj x in Eq. ([2]). Typical experimental parameters range from 
qj x = 2n x 10 Hz to 2tt x 75 Hz fl-i]. 

To get a impression of the energy scales, we show some val- 
ues for the Hubbard parameters obtained by our band structure 
approach in Tab. [IV] By setting si = the on-site energies are 
solely due to the additional harmonic potential. At the outer 
rims of the lattice (sites 1 and 30) they have the value 6 max . 

Up to (jj x - 2n x 25 Hz, e max /J is an order of magnitude 
smaller than U/J. For this reason, the phase diagram remains 
practically unaltered between co x = Hz and 2n x 25 Hz as 
can be seen by comparing Figs. [Sa)-(cl) and I3a)-(d). For 
(jj x - 2n x 50 Hz Tab. HV1 shows e max /J is still about a fac- 
tor 3 smaller than U/J. As a consequence the onset of the 
BG phase in Figs. [T^e) -(h) already appears at s\ « 0.4 instead 
of s\ « 0.6 for oj x = Hz. Besides the earlier onset of the 
BG phase also its gross structure changes. The lobe around 
s\ = 0.6 in Fig. [7{e) is suppressed compared to the calcula- 
tions for co x < 2n x 50 Hz . Also the maximum fluctuations 
indicate that the redistribution of particles becomes smoother. 
This is because for co x = Hz the superlattice topology ex- 
hibits only 5 different on-site energies. With the additional 
harmonic potential the number of different on-site energies in- 
creases and, therefore, the extended domains in the BG phase 
shrink. For co x = 2n x 75 Hz the parameters e m3iX /J and U/J 
become comparable and the phase diagram changes dramati- 
cally. In Fig. |7ti) the homogeneous MI domain shrinks to a 
small region (s2 = 12-16 and s\ = - 0.2). Furthermore, 
a clear detection of the BG phase becomes difficult since the 
characteristic increase of the visibility along the MI to BG 
transition is no longer visible in Fig.^l). TableHVlreveals that 
at cj x = 2nx 100 Hz the on-site energies 6 max clearly dominate 
the energy scale. Thus, even the transition from the SF phase 
to the MI phase is no longer observable in the investigated 
parameter range of S2. The MI domain has completely disap- 
peared in Fig. [3m) and the visibility remains large throughout 
the whole range of ^2 for s\ < 0.2. 
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FIG. 7: (color online) Phase diagram in terms of energy gap AE/J, condensate fraction f c , maximum number fluctuation <x max , and visibility 
y for I = N = 30, a s = 109 r Bo hr, w ± - 2n x 17 kHz. First row: oj x - 2n X 25 Hz, second row: a> x - 2n X 50 Hz, third row: a> x = 2nx 75 Hz, 
fourth row: ex) x = In X 100 Hz. 



From the above discussion we conclude that the smaller the 
longitudinal trapping frequency co x , the easier is a clear dis- 
tinction between the SF, MI and BG phases. Thus, any exper- 
iment with a focus on the phase diagram of ultracold atoms in 
an optical superlattice should be designed such that the longi- 
tudinal trapping frequency is kept small. 



VIII. TRANSVERSE TRAPPING FREQUENCY cj ± 

Finally we study the dependence of the (s\, s^) phase dia- 
gram on the intensity of the transverse lasers through a varia- 
tion of the transverse trapping frequency oj ± . For the sake of 
simplicity we assume co x = Hz. 

From Eq. © it follows that the interaction energy Ut is 
proportional to oj ± while 6/ and J\ are independent of a) ± . A 
larger value of oj ± will, therefore, shift the SF to MI transition 
towards smaller S2. Also the spread of the on-site energies 
must increase to overcome energy cost of a double occupancy 
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FIG. 8: (color online) Same parameter set as in Fig.[6]but with transverse trapping frequency a> ± = 2n x 40 kHz. While the gross structure of 
the phase diagram is independent of oj ± , the scales of the s 2 and s\ axes change. 



and consequently the MI to quasi BG transition will shift to- 
wards larger s\. In Fig. [8] we show the phase diagrams for 
oj ± = 2n x 40 kHz and all other parameters unchanged. The 
gross structure of the phase diagram remains the same. How- 
ever, in accordance with our considerations, the energy gap 
is AE/J = 60 already at ^2 = 12.5 instead of ^2 = 16 for 
o) ± = 2n x 17 kHz (both s\ = 0). Furthermore, the transi- 
tion from the homogeneous MI to the quasi BG phase occurs 
around s\ « 1.5 as compared to S2 = 0.6 for oj ± = 2n x 17 
kHz. 

As a consequence of the dependence of the (s^, s\) phase di- 
agram on the transverse trapping frequency, the precise iden- 
tification of the phase boundaries is intimately connected to a 
well defined value of a) ± . 



IX. SUMMARY & CONCLUSIONS 

We have studied the experiment- specific phase diagram of 
ultracold 87 Rb atoms in an one-dimensional two-color super- 
lattice with respect to the parameters of the experiment. Band 
structure calculations were employed to obtain the generic pa- 
rameters of the Hubbard model from the experiment- specific 
parameters. These band structure calculations were also used 
to confirm the applicability of the Hubbard model in the in- 
vestigated parameter range. 

In order to solve the many-body problem for realistic lat- 
tice lengths and particle numbers we have used the density- 
matrix renormalization-group algorithm. Through a thorough 



benchmark of our DMRG calculations we demonstrated that 
all observables are perfectly converged and can practically be 
considered as exact solutions of the many-body problem. Fur- 
thermore, a detailed finite-size analysis for all observables has 
underlined the significance of our results for realistic experi- 
mental system sizes. 

Our calculations of the phase diagrams show that all rel- 
evant quantum phases can be addressed by only varying the 
intensities of the two lasers that generate the optical superlat- 
tice. For a longitudinal trapping frequency co x < 25 Hz all dif- 
ferent quantum phases can be clearly distinguished by means 
of the presented observables. However, larger values of the 
longitudinal trapping frequency lead to radical changes in the 
structure of the phase diagram and make a clear identification 
of the Bose-glass phase impossible. 

We also showed that the gross structure of the phase dia- 
gram does not depend on the transverse trapping frequency 
o) ± , i.e., the intensity of the lasers in the directions perpendic- 
ular to the ID lattice. However, due to the linear dependence 
of the interaction energy on the transverse trapping frequency, 
the position of the transition lines in the phase diagram cru- 
cially depend on that parameter. 
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